! Bessel I_0(x) function in double precision
!
      function dbesi0(x)
      implicit real*8 (a - h, o - z)
      dimension a(0 : 64), b(0 : 69), c(0 : 44)
      data (a(i), i = 0, 12) / 
     &    8.5246820682016865877d-11, 2.5966600546497407288d-9, 
     &    7.9689994568640180274d-8, 1.9906710409667748239d-6, 
     &    4.0312469446528002532d-5, 6.4499871606224265421d-4, 
     &    7.9012345761930579108d-3, 7.1111111109207045212d-2, 
     &    4.4444444444472490900d-1, 1.7777777777777532045d0, 
     &    4.0000000000000011182d0, 3.9999999999999999800d0, 
     &    1.0000000000000000001d0 / 
      data (a(i), i = 13, 25) / 
     &    1.1520919130377195927d-10, 2.2287613013610985225d-9, 
     &    8.1903951930694585113d-8, 1.9821560631611544984d-6, 
     &    4.0335461940910133184d-5, 6.4495330974432203401d-4, 
     &    7.9013012611467520626d-3, 7.1111038160875566622d-2, 
     &    4.4444450319062699316d-1, 1.7777777439146450067d0, 
     &    4.0000000132337935071d0, 3.9999999968569015366d0, 
     &    1.0000000003426703174d0 / 
      data (a(i), i = 26, 38) / 
     &    1.5476870780515238488d-10, 1.2685004214732975355d-9, 
     &    9.2776861851114223267d-8, 1.9063070109379044378d-6, 
     &    4.0698004389917945832d-5, 6.4370447244298070713d-4, 
     &    7.9044749458444976958d-3, 7.1105052411749363882d-2, 
     &    4.4445280640924755082d-1, 1.7777694934432109713d0, 
     &    4.0000055808824003386d0, 3.9999977081165740932d0, 
     &    1.0000004333949319118d0 / 
      data (a(i), i = 39, 51) / 
     &    2.0675200625006793075d-10, -6.1689554705125681442d-10, 
     &    1.2436765915401571654d-7, 1.5830429403520613423d-6, 
     &    4.2947227560776583326d-5, 6.3249861665073441312d-4, 
     &    7.9454472840953930811d-3, 7.0994327785661860575d-2, 
     &    4.4467219586283000332d-1, 1.7774588182255374745d0, 
     &    4.0003038986252717972d0, 3.9998233869142057195d0, 
     &    1.0000472932961288324d0 / 
      data (a(i), i = 52, 64) / 
     &    2.7475684794982708655d-10, -3.8991472076521332023d-9, 
     &    1.9730170483976049388d-7, 5.9651531561967674521d-7, 
     &    5.1992971474748995357d-5, 5.7327338675433770752d-4, 
     &    8.2293143836530412024d-3, 6.9990934858728039037d-2, 
     &    4.4726764292723985087d-1, 1.7726685170014087784d0, 
     &    4.0062907863712704432d0, 3.9952750700487845355d0, 
     &    1.0016354346654179322d0 / 
      data (b(i), i = 0, 13) / 
     &    6.7852367144945531383d-8, 4.6266061382821826854d-7, 
     &    6.9703135812354071774d-6, 7.6637663462953234134d-5, 
     &    7.9113515222612691636d-4, 7.3401204731103808981d-3, 
     &    6.0677114958668837046d-2, 4.3994941411651569622d-1, 
     &    2.7420017097661750609d0, 14.289661921740860534d0, 
     &    59.820609640320710779d0, 188.78998681199150629d0, 
     &    399.87313678256011180d0, 427.56411572180478514d0 / 
      data (b(i), i = 14, 27) / 
     &    1.8042097874891098754d-7, 1.2277164312044637357d-6, 
     &    1.8484393221474274861d-5, 2.0293995900091309208d-4, 
     &    2.0918539850246207459d-3, 1.9375315654033949297d-2, 
     &    1.5985869016767185908d-1, 1.1565260527420641724d0, 
     &    7.1896341224206072113d0, 37.354773811947484532d0, 
     &    155.80993164266268457d0, 489.52113711585409180d0, 
     &    1030.9147225169564806d0, 1093.5883545113746958d0 / 
      data (b(i), i = 28, 41) / 
     &    4.8017305613187493564d-7, 3.2613178439123800740d-6, 
     &    4.9073137508166159639d-5, 5.3806506676487583755d-4, 
     &    5.5387918291051866561d-3, 5.1223717488786549025d-2, 
     &    4.2190298621367914765d-1, 3.0463625987357355872d0, 
     &    18.895299447327733204d0, 97.915189029455461554d0, 
     &    407.13940115493494659d0, 1274.3088990480582632d0, 
     &    2670.9883037012547506d0, 2815.7166284662544712d0 / 
      data (b(i), i = 42, 55) / 
     &    1.2789926338424623394d-6, 8.6718263067604918916d-6, 
     &    1.3041508821299929489d-4, 1.4282247373727478920d-3, 
     &    1.4684070635768789378d-2, 1.3561403190404185755d-1, 
     &    1.1152592585977393953d0, 8.0387088559465389038d0, 
     &    49.761318895895479206d0, 257.26842323135291380d0, 
     &    1066.8543146269566231d0, 3328.3874581009636362d0, 
     &    6948.8586598121634874d0, 7288.4893398212481055d0 / 
      data (b(i), i = 56, 69) / 
     &    3.4093503681970328930d-6, 2.3079025203103376076d-5, 
     &    3.4691373283901830239d-4, 3.7949949772229085450d-3, 
     &    3.8974209677945602145d-2, 3.5949483804148783710d-1, 
     &    2.9522878893539528226d0, 21.246564609514287056d0, 
     &    131.28727387146173141d0, 677.38107093296675421d0, 
     &    2802.3724744545046518d0, 8718.5731420798254081d0, 
     &    18141.348781638832286d0, 18948.925349296308859d0 / 
      data (c(i), i = 0, 8) / 
     &    2.5568678676452702768d-15, 3.0393953792305924324d-14, 
     &    6.3343751991094840009d-13, 1.5041298011833009649d-11, 
     &    4.4569436918556541414d-10, 1.7463930514271679510d-8, 
     &    1.0059224011079852317d-6, 1.0729838945088577089d-4, 
     &    5.1503226936425277380d-2 / 
      data (c(i), i = 9, 17) / 
     &    5.2527963991711562216d-15, 7.2021184814210056410d-15, 
     &    7.2561421229904797156d-13, 1.4823121466731042510d-11, 
     &    4.4602670450376245434d-10, 1.7463600061788679671d-8, 
     &    1.0059226091322347560d-6, 1.0729838937545111487d-4, 
     &    5.1503226936437300716d-2 / 
      data (c(i), i = 18, 26) / 
     &    1.3365917359358069908d-14, -1.2932643065888544835d-13, 
     &    1.7450199447905602915d-12, 1.0419051209056979788d-11, 
     &    4.5804788198059832600d-10, 1.7442405450073548966d-8, 
     &    1.0059461453281292278d-6, 1.0729837434500161228d-4, 
     &    5.1503226940658446941d-2 / 
      data (c(i), i = 27, 35) / 
     &    5.3771611477352308649d-14, -1.1396193006413731702d-12, 
     &    1.2858641335221653409d-11, -5.9802086004570057703d-11, 
     &    7.3666894305929510222d-10, 1.6731837150730356448d-8, 
     &    1.0070831435812128922d-6, 1.0729733111203704813d-4, 
     &    5.1503227360726294675d-2 / 
      data (c(i), i = 36, 44) / 
     &    3.7819492084858931093d-14, -4.8600496888588034879d-13, 
     &    1.6898350504817224909d-12, 4.5884624327524255865d-11, 
     &    1.2521615963377513729d-10, 1.8959658437754727957d-8, 
     &    1.0020716710561353622d-6, 1.0730371198569275590d-4, 
     &    5.1503223833002307750d-2 / 
      w = abs(x)
      if (w .lt. 8.5d0) then
          t = w * w * 0.0625d0
          k = 13 * (int(t))
          y = (((((((((((a(k) * t + a(k + 1)) * t + 
     &        a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + 
     &        a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + 
     &        a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + 
     &        a(k + 11)) * t + a(k + 12)
      else if (w .lt. 12.5d0) then
          k = int(w)
          t = w - k
          k = 14 * (k - 8)
          y = ((((((((((((b(k) * t + b(k + 1)) * t + 
     &        b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + 
     &        b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + 
     &        b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + 
     &        b(k + 11)) * t + b(k + 12)) * t + b(k + 13)
      else
          t = 60 / w
          k = 9 * (int(t))
          y = ((((((((c(k) * t + c(k + 1)) * t + 
     &        c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + 
     &        c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + 
     &        c(k + 8)) * sqrt(t) * exp(w)
      end if
      dbesi0 = y
      end
!
